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We introduce numerical algorithms for initializing multidimensional simulations of stellar explosions with ID stellar evolution models. 
The initial mapping from ID profiles onto multidimensional grids can generate severe numerical artifacts, one of the most severe of which 
is the violation of conservation laws for physical quantities. We introduce a numerical scheme for mapping ID spherically-symmetric 
data onto multidimensional meshes so that these physical quantities are conserved. We validate our scheme by porting a realistic ID 
Lagrangian stellar profile to the new multidimensional Eulerian hydro code CASTRO. Our results show that all important features in the 
profiles are reproduced on the new grid and that conservation laws are enforced at all resolutions after mapping. We also introduce a 
numerical scheme for initializing multidimensional supernova simulations with realistic perturbations predicted by ID stellar evolution 
models. Instead of seeding 3D stellar profiles with random perturbations, we imprint them with velocity perturbations that reproduce the 
Kolmogorov energy spectrum expected for highly turbulent convective regions in stars. Our models return Kolmogorov energy spectra 
and vortex structures like those in turbulent flows before the modes become nonlinear. Finally, we describe approaches to determining the 
resolution for simulations required to capture fluid instabilities and nuclear burning. Our algorithms are applicable to multidimensional 
simulations besides stellar explosions that range from astrophysics to cosmology. 
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1 Introduction 

Multidimensional simulations shed light on how fluid instabilities arising in supern ovae (SNe) mix ejecta 
dHerant and Wooslevl Il994l . I.Toggerst et all 120091 . l2010l . iJoggerst and WhaleiilEoilh . Unfortunately, com- 
puting fully self-consistent 3D stellar evolution models, from their formation to collapse, for the explosion 
setup is still beyond the realm of contemporary computational power. One alternative is to hrst evolve the 
main sequence star in a ID stellar evolution code in which th e equations of mom entum, ene rgy and mass 
are s olved on a spherically-symmetric grid, such as KEPLER ( Weaver et al. 19781 ) or MESA ( Paxton et al. 
201 lh . Once the star reaches the pre-su pernova phase, its ID profiles can then be mapped into multidi 



mens ional hydro codes such as CASTRO ( Almgren et aZ.I [201ol . Izhang et all 1201 ll ) or FLASH ( Fryxell et al 
2000) and continue to be evolved until the star explodes 



Differences between codes in dimensionality and coordinate mesh can lead to numerical issues such as 
violation of conservation of mass and energy when profiles are mapped from one code to another. A first, 
simple approach could be to initialize multidimensional grids by linear interpolation from corresponding 
mesh points on the ID profiles. However, linear interpolation becomes invalid when the new grid fails 
to resolve critical features in the original profile such as the inner core of a star. This is especially true 
when porting profiles from ID Lagrangian codes, which can easily resolve very small spatial features in 
mass coordinate, to a fixed or adaptive Eulerian grid. Besides conservation laws, some physical processes 
such as nuclear burning are very sensitive to temperature, so errors in mapping can lead to very different 
outcomes for the simulations such as altering the nucleosynt hesis and energetics of SNe. Only few studies 



have examined mapping ID profiles to 2D or 3D meshes ( Zingale et a l. 2002), and none address the 



conservation of physical quantities by such procedures. We examine these issues and introduce a new 
scheme for mapping ID data sets to multidimensional grids. 
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Seeding the pre-supernova profile of the star with realistic perturbations is important to understanding 
how fluid instabilities later erupt and mix the star during the e xplosion. Massive stars u sually develop 
convective zones prior to exploding as SNe (jWooslev et of J [2002. Heger and Wooslevll2002h . Multidimen- 
sional stellar evolution models suggest that the fluid ins ide the convective regions can be highly turbulent 
( Porter and Woodward! 2000l . Arnett and Meakin 201 ll ). However, in lieu of the 3D stellar evolution cal- 
culations necessary to produce such perturbations from first principles, multidimensional simulations are 
usually just seeded with random p erturbations. In reality, if the star is convective and the fluid in those 
zones is turbulent (|Davidsonl 12004 ) . a better approach is t o imprint th e multidimensional profiles with 
velocity perturbations with a Kolmogorov energy spectrum ( Frisch 19951 ). 

Besides implementing realistic initial conditions, care must be taken to determine the resolution multi- 
dimensional simulations require to resolve the most important physical scales and yield consistent results 
given the computational resources that are available. We provide a systematic approach for finding this 
resolution for multidimensional stellar explosions. The structure of paper is as follows. In § [2] we describe 
the key features of the KEPLER and CASTRO codes. We describe our initial mapping scheme and demonstrate 
it by porting a massive star model from KEPLER to CASTRO in § We review our scheme for seeding 2D and 
3D stellar profiles with turbulent perturbations and present hydrodynamic simulations done with these 
profiles in in CASTRO in § UJ We provide a strategy for finding the proper resolution for multidimensional 
simulations with multiscale processes such as hydrodynamics and nuclear burning in § [5] and conclude the 
results in § [6j 



2 Stellar Model 



We model the evolution of main sequence stars with KEPLER ( Weaver et al. 19781 ). a ID Lagrangian stellar 
evolution code. KEPLER solves the evolution equations for mass, momentum, and energy, including relevant 
physical processes such as nuclear burning and mixing due to convection. When the star reaches the pre- 
supernova phase (hundreds of seconds prior to launching the SN shock), we map its ID profiles onto a 
multidimensional grid in CASTRO. When the star explodes, its initial spherical symmetry is broken by fluid 
instabilities formed during the explosion that cannot be modeled by ID calculations. Hence, we follow the 
evolution of the star in CASTRO until it explodes. 

CASTRO ( Almgren et ai]|201Ct IZhang et aZ.ll201lh is a massively parallel, multidimensional Eulerian adap- 
tive mesh refinement (AMR) radiation-hydrodynamics code for astrophysical applications. Its time in- 
tegration of the hydrodynamics equations is based on a higher-order, unsplit Godunov scheme. Block- 
structured AMR with subcycling in time applies high spatial resolution to where it is needed most. We 
use the Helmholtz equation of state ( Timmes and Swestv 2000l ) with density, temperature, and elemental 
abundances; it includes contributions by non-degenerate and degenerate relativistic and non-relativistic 
electrons, electron-positron pair production, ions and radiation. The gravitational field is calculated with 
a monopole approximation derived from a radial average of the den sity on the multidimensional grid . 
We have implemented several reaction networks (7, 13, 19 isotopes) ( Weaver et al. 19781 . Timmesl 1999) 
in CASTRO for calculating nucleosynthesis and energetics in thermonuclear SNe. The most comprehensive 
network includes alpha-chain reactions, heavy-ion reactions, hydrogen burning cycles, photo-disintegration 
of heavy elements, and energy loss by neutrinos. 



3 Conservative Mapping 

Since the star is very nearly in hydrostatic equilibrium and we want to conserve total energy, care must be 
taken when mapping its profile from the uniform Lagrangian grid in mass coordinate to the new Eulerian 
spatial grid. Our method preserves the conservation of quantities such as mass and energy on the new 
mesh that are analytically conserved in the evolution equations. Although this reconstruction does not 
guarantee that the star will be hydrostatic, it is a physically motivated constraint and sufficient for our 
simulations. The algorithm we describe is specific to our models but can be easily generalized to mappings 
of other ID data to higher dimensional grids. 
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Figure 1. Constructing a piecewise-linear conservative profile: The rectangular bins represent the original ID profile. The areas of 
different colors represent conserved quantities such as mass and internal energy. The conservative profile connects adjacent bins and 
makes a =a' = ^Hx min(A,B), c = c', and b = b'. Note that uniform zones in mass (Lagrangian coordinate) lead to nonuniform bins 
in volume coordinate, as shown above. 



3.1 Method 



First, we construct a continuous (C°) function that conserves the physical quantity upon mapping onto 
the new grid. An ideal choice for interpolation is the volume coordinate V, the volume enclosed by a given 
radius from the center of the star. Then, integrating a density px (which can represent mass or internal 
energy density) with respect to the volume coordinate yields a conserved quantity X 



X 



v 2 



Px dV, 



(1) 



such as the total mass or total internal energy lying in the shell between V% and V%- 

Next, we define a piecewise linear function in volume V that represents the conserved quantity px, 
preserves its monotonicity (no new artificial extrema) , and is bounded by the extrema of the original data. 
The segments are constructed in two stages. First, we extend a line across the interface between adjacent 
zones that either ends or begins at the center of the smaller of the two zones, as shown in Figure Q] (note 
that uniform zones in mass coordinate do not result in uniform zones in V). The slope of the segment is 
chosen such that the area trimmed from one zone by the segment (a and b) is equal to the area added 
under the segment in the neighboring bin (a'=a and b'=b). 

If the two segments bounding a and a', and b and b' are joined together by a third in the center zone 
in Figure [H two "kinks", or changes in slope, can arise in the interpolated quantity there; plus, the slope 
of the fiat central segment is usually a poor approximation to the average gradient in that interval. We 
therefore construct two new segments that span the entire central zone and connect with the two original 
segments where they cross its interfaces, as shown in Figure [U The new segments join each other at the 
position in the central bin where the areas c and c' enclosed by the two segments are equal (note that they 
in general have different slopes). After repeating this procedure everywhere on the grid, each bin will be 
spanned by two linear segments that represent the interpolated quantity px at any V within the bin and 
have no more than one kink in px across the zone. Our scheme introduces some smearing (or smoothing) 
of the data, but it is limited to at most the width of one zone on the original grid. 

The result of our interpolation scheme is a piecewise linear reconstruction in V of the original profile in 
mass coordinate for which the quantity px can be determined at any V, not just the radii associated with 
the zone boundaries in the Lagrangian grid. We show this profile as a function of the radius associated 
with the volume coordinate V for a zero-metallicity 200 M star with r ~ 2 x 10 13 cm from KEPLER 
dHeger and Wooslevl I2OO2L l20ld ). 

We populate the new multidimensional grid with conserved quantities from the reconstructed stellar 
profiles as follows. First, the distance of the selected mesh point from the center of the new grid is 
calculated. We then use this radius to obtain its V to reference the corresponding density in the piecewise 
linear profile of the star. The density assigned to the zone is then determined from adaptive iterative 
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Mc=Vc*p c 



Md=Vd-pd 




Mi = Ma+Mb+Mc+...+Mh 



Figure 2. Volume subsampling: We first use the center of volume element (Vo) to obtain its density po from the conservative profile and 
then calculate its mass Mo = Vo x po . We then partition the original volume element as shown and calculate the mass of each subvolumc 
in the same manner as Mo- We obtain Mi by summing over eight subvolumes M a ,Mj,,M c ,...,M^. We then compare Mi and Mo; if their 
relative error is greater than some predetermined tolerance the process is recursively repeated until |(Mj-Mj_i)/Mj| is less then 10 — 4 . 

subsampling. This is done by first computing the total mass of the zone by multiplying its volume by the 
interpolated density. We then divide the zone into equal subvolumes whose sides are half the length of 
the original zone. New V are computed for the radii to the center of each of these subvolumes and their 
densities are again read in from the reconstructed profile. The mass of each subvolume is then calculated by 
multiplying its interpolated density by its volume element (see Figure [2]). These masses are then summed 
and compared to the mass previously calculated for the entire cell. If the relative error between the two 
masses is larger than the desired tolerance, each subvolume is again divided as before, masses are computed 
for all the constituents comprising the original zone, and they are then summed and compared to the zone 
mass from the previous iteration. This process continues recursively until the relative error in mass between 
the two most recent consecutive iterations falls within an acceptable value, typically 10 -4 . The density we 
assign to the zone is just this converged mass divided by the volume of the entire cell. This method is used 
to map internal energy density and the partial densities of the chemical species to every zone on the new 
grid. The total density is then obtained from the sum of the partial densities; pressure and temperature 
in turn are determined from the equation of state. This method is easily applied to hierarchy geometry of 
the target grid. 



We port a ID stellar model from KEPLER into CASTRO to verify that our mapping is conservative. As an 
example, we use a 200 M zero-metalicity pre-supernova star. 

As we show in Figure O our piecewise linear fits to the KEPLER data reproduce the original stellar profile. 
Because our fits smoothly interpolate the block histogram structure of the KEPLER bins (especially at larger 
radii), they reduce the number of unphysical sound waves that would have been introduced in CASTRO by 
the discontinuous interfaces between these bins in the original dats0- The density profile is key to the 
hydrodynamic and gravitational evolution of the explosion, and the temperature profile is crucial to the 
nuclear burning that powers the explosion. 

We first map the profile onto a ID grid in CASTRO and plot the mass of the star as a function of grid 
resolution in Figure HI The mass is independent of resolution for our conservative mapping because we 
subsample the quantity in each cell prior to initializing it, as described above. In contrast, the total mass 
from linear interpolation is very sensitive to the number of grid points but does eventually converge when 
the number of zones is sufficient to resolve the core of the star, in which most of its mass resides. 

We next map the KEPLER profile onto a 2D cylindrical grid (r, z) and a 3D cartesian grid (x, y, z) in 
CASTRO. The only difference between mapping to ID, 2D, and 3D is the form of the volume elements used 



3.2 Results 



*1D data usually provides zone-averaged values, hence a continuous and conservative profile needs to be reconstructed from zone-averaged 
values. 
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(b) Inner temperature profile 



Figure 3. Inner density and temperature profiles of a 200 Mq presupernova: Our piecewise linear profiles (green lines) fit their original 
KEPLER model (red crosses) very well. Since we map internal energy (a conserved quantity) rather than temperature, we calculate T from 
the equation of state using the density, element abundance, and internal energy. 
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Figure 4. Total mass of the star on the ID CASTRO grid as a function of number of zones: Conservative mapping (blue) preserves the 
mass of the star at all resolutions, while linear interpolation (orange) converges to 200 Mq at a resolution of ~ a few xlO 4 , when the 
grid begins to resolve the core of the star (~ 10 9 cm). Even at very high resolutions, the results of linear interpolation are still off by a 
few percent from the targeted mass and start to be saturated at ~ 10 5 zones because the linear interpolation profile is a non-conservative 
profile in nature. 



to subsample each cell, which are Airr 2 dr, 2Trrdrdz, dxdydz, respectively. We show the mass of the star as a 
function of resolution in Figure 5(a) Conservative mapping again preserves its mass at all grid resolutions. 



In 2D, more zones are required for linear interpolation to converge to the mass of the star. To further 
validate our conservative scheme, we map just the helium core of the star (~ 100 Mq with r ~ 10 10 cm) 
onto the 2D grid. The helium core is crucial to modeling thermonuclear supernovae because it is where 
explosive burning begins. We show its mass as a function of resolution in Figure 5(b) We again recover 
all the mass of the core at all resolutions whereas linear interpolation overestimates the mass by at least 
~ 1 %, even with large numbers of zones. 

Because of the property of reconstruction, conservative mapping is effective in 3D but requires much 
more computational time to subsample each cell to convergence. Furthermore, an impractical number 
of zones is needed for linear interpolation to reproduce the original mass of the star. 3D results is not 
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Number of Zones in R and Z Number of Zones in R and Z 

(a) 2D Mapped Mass of Entire Star (b) 2D Mapped Mass of Helium Core 

Figure 5. (a) Total mass of the star on the new 2D CASTRO grid as a function of number of zones in both r and z: Conservative mapping 
(blue) recovers the mass of the star at all resolutions and linear interpolation (orange) approaches 200 Mq at a resolution of ~ 2048 2 . 
(b) Total mass of the He core on the 2D CASTRO grid as a function of number of zones in both r and z: conservative mapping (blue) 
preserves its original mass at all resolutions while linear interpolation (orange) begins to converge to 100 Mq at a resolution of 64 2 , but 
it is still off by ~ 1% even as the resolution approaches ~ 2048 2 because the linear interpolation profile is a non-conservative profile in 
nature. 




Figure 6. Convective core of star: The ID stellar model with a stellar radius, R can predict the size of inner convective core, r shown as 
the wavy meshes above as well as the fluid velocities inside it. 

presented here. We note that our method also works with AMR grids because both V and the interpolated 
quantities can be determined, and subsampling can be performed on every grid in the hierarchy. For the 
given domain, the results of conservative mapping are independent of the levels of AMR. 

4 Initial Perturbation 

Next we describe our scheme for seeding 2D and 3D stellar profiles with turbulent perturbations and 
present stellar evolution simulations with CASTRO with these profiles. In our setup, the perturbations have 
following properties (see Figure [6]) : 

(i) The perturbations are imprinted in the gas velocity, and their net momentum flux must be zero. 

(ii) They are seeded in convectively unstable regions with a velocity spectrum V(k) ~ /c _5//6 , where k 
is the wave number and the power index —5/6 is for a Kolmogorov spectrum with an assumption 
of constant density. 
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4.1 2D Perturbation 

We first consider the mapping onto a polar coordinate grid in r and 9. To enforce zero net momentum and 
the boundary conditions in the simulation, we define a new variable x = 1 + cos 9 instead of using 9. The 
momentum flux of a density p and velocity v r is then 

/•7T 1-2 

I 2nr 2 pv r sin 9d9 = 2vrr 2 pv r dx = (2) 
Jo Jo 

if v r has the form cos(27rnx), where n is an integer. When 9 = 0, ir (the boundaries of a 2D grid), x = 
2, yields the maximum values for v r that satisfy the boundary conditions in 2D cylindrical coordinates 
in CASTRO. There are two ph ysical scales that co nstrain the wavelength of the perturbation in r. Based 
on the mixing length theory (|Cox and Giulil l 19681 ). the eddy size of turbulence is a x H p , a is the mixing 



length parameter, H p is the pressure scale height. Here, we set a = 1.0. Since the perturbation only seeded 
in the convective zones, it is confined, inside domain D c — v u - rb, where r u and are its upper and 
lower boundaries. The maximum wavelength of the perturbation must be smaller than D c and H p . Inside 
a convective zone, we define a new variable, y = j/^r) • ^ e a ^ so define two oscillatory functions in x 
and y to generate the circular patterns that mimic the vortices of a turbulent fluid. Since the fluid inside 
the convective zone is turbulent, its energy spectrum E(k) ~ A; -5 / 3 . Assuming a constant density, the 
corresponding velocity spectrum is V{k) ~ /c~ 5//6 . The perturbed velocity then has the form, 

V per b t r(x, y) = - Y, V p ■ cos(27rax) • cos(27r6y + a b ), 

a b 

V per b,e(x, y) = Y, Y, V p ' sin(27rax) • cos(2irby + a b ), (3) 
V p = V d (r)b- 5 /i 



where V per b, r an d V pe rbfi are th e perturbed velocities in the r and 9 directions and a and b are angular and 
radial wavenumbers. ID models only provide the information of convective velocities, Vd(r) along radial 
direction which can be treated as average velocities of angular directions, so we scale the amplitude of 
the perturbed velocity based on radial wavenumber b. We also use a random phase, ab, to smooth out 
numerical discontinuities caused by the perturbed modes while summing. Equations ([3]) by construction 
satisfy V per bfi(r, 9) = when 9 = and it, the boundary conditions in 9 on the 2D grid. Besides, we assume 
there is no overshoting between the upper and lower boundaries of convective zones and enforce V per b t r to 
zero by setting ab = tt/2 when is an integer. The ultimate wavenumbers of a and b are also limited 
by D c , H p , and the resolution of simulation, H res . 



4.2 3D Perturbation 

In 3D, we use spherical coordinate r, 9, and (ft. Similar to 2D, we construct an oscillatory function for (9, 4>) 
by using spherical harmonics, Y^ m (9, (ft), where I and m are the wavenumbers in 9 and (j). If the velocities in 
the form of Y[ jm (9,(f)), they automatically conserve momentum flux while summing all the modes l,m. In 
the radial direction, we use cos(cy), where c is the wavenumber in the radial direction and y is as defined 
in 2D. The perturbation then has the following form, 



V per b,x{r, 0,4>) = Vp er b sin(6>) cos(^), 
Vperb,y(r,9,<j)) = V perb sin(0) sin(^) 5 

V per b, z {r, 9, <{>) = V perb cos(9), ( 4 ) 

Vperb = EEE l p- Yl,m{9 + ^Ira, 4> + <^lm) ' COs(27TCy + A c ), V p = V d (r)c~ 5/6 , 
c I rn 

where V per b,x, Vperb,y, Vperb,z are the perturbed velocities in the x, y, and z directions. We sum over the 
modes, applying random phases w/ m and A c to smooth out numerical discontinuities caused by different 
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perturbed modes. Similar to 2D, V p is only scaled by radial wavenumber c. Because there is no reflective 
boundary conditions for 3D, we only take care of the boundary conditions in radial direction, we again 
assume there is no overshoting outside the boundaries of convective zones, so enforce V per b to zero by 
setting A c = ir/2 when "cy" is an integer. 



4.3 Results 



We first initialize perturbations on a 2D grid with a profile that is derived from a ID KEPLER stellar evolutio n 
calculation. The perturbations are confined to regions that are convectively unstable ( Heger et al. 2000l ). 



For the magnitude V^(r) of the perturbed velocity adopts the diffusion velocity, which is usually ~ 1 - 10 % 
of the local sound speed. We again consider a zero-metalicity 200 M star in the pre-supernova phase. This 
star develops a large convection zone that extends out to the hydrogen envelope. We show the magnitude 
of the perturbed velocity generated by the two oscillatory functions discussed above on our 2D grid in 



Figure 7(a) The velocity field satisfies the reflecting boundary conditions on the 2D grid at = and 7r. 
In the right panel we show velocity vectors in the selected subregion on the left (blue rectangle). A clear 
vortex pattern that mimics a turbulent fluid is clearly visible. Figure |7(b)| shows the energy spectrum of 
the fluid, which is basically a Kolmogorov spectrum E(k) ~ /c -5 / 3 except for fluctuations in part caused 
by the random phases in the sum over modes in r and Vd(r) is not a constant across the convective region 
that produces an offset in the smaller k region. The energies would converge to the Kolmogorov spectrum 
in the limit of large k, but the maximum k of our simulation is limited to the resolution of the grid. 



We next port our ID KEPLER model to a 3D grid. In Figure 8(a) we show a slice of the magnitude of the 
perturbed velocity, which again exhibits the clear cell pattern reminiscent of the vortices of a turbulent 
fluid. The velocity pattern in 3D is more irregular than in 2D. We show the energy spectrum of the velocity 
field in Figure [8(b)[ which is similar to that of our 2D spectrum but with larger fluctuations that are again 
due to the random phases we assign to each spherical harmonic and and the Vd(r) is not a constant 
across the convective region that produces an offset in the smaller k region. Our 2D and 3D examples 
demonstrate that our scheme effectively generates turbulence fluid perturbations like those found in the 
convective regions of massive stars, with the desired initial velocity patterns and energy power spectra. 

We do not claim the models here can fully reproduce the true turbulence found in simulations or labo- 
ratories. The scheme here is the first attempt to model the initial perturbations based on a more realistic 
setup. Seeding initial perturbation is to trigger the fluid instabilities on multidimensional simulations , so 
we can study how they evolve with their surroundings. When the fluid instabilities start to evolve nonlin- 
early, the initial imprint of perturbation would be smeared out. The random perturbations and turbulent 
perturbations then give consistent results. Depending on the problems, the random perturbation might 
take a longer time to evolve the fluid instabilities into turbulence because more relaxation time is required. 



5 Resolving the Early Stages of the Explosion 

Simulations that include nuclear burning, which governs nucleosynthesis and the energetics of explosion, 
are very different from purely hydro dynamical models because of the more stringent resolution required 
to resolve the scales of nuclear burning and the onset of fluid instabilities in the simulations. Because 
energy generation rates due to burning are very sensitive to temperature, errors in these rates as well as in 
nucleosynthesis can arise in zones that are not fully resolved. We determine the optimal resolution with a 
grid of ID models in CASTRO. Beginning with a crude resolution, we evolve the pre-supernova star and its 
explosion until all burning is complete and then calculate the total energy of the supernova, which is the 
sum of the gravitational energy, internal energy, and kinetic energy. We then repeat the calculation with 
the same setup but a finer resolution and again calculate the total energy of the explosion. We repeat this 
process until the total energy is converged. As shown in Figure [9j our example of 200 M presupernova 
converges when the resolution of the grid approaches 10 8 cm. 

The time scales of burning (dtb) and hydro dynamical (dth) can be very disparate, so we adopt time 
steps of min(dth, dtf,) in our simulations, where dth = d f\ v \ ; dx is the grid resolution, c s is the local sound 
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Figure 7. (a) 2D perturbed velocity in the interior of the star on physical scales ~ 10 12 cm: the closeup is the velocity vector field 
corresponding to the blue rectangle and exhibits a vortex pattern similar to that of a turbulent fluid, (b) Normalized kinetic energy 
power spectrum of a 2D perturbed field: The dotted red line is the Kolmogorov spectrum, E(k) ~ fc -5 / 3 . The peak of the Kolmogorov 
spectrum is adjusted to fit the data. The scale of H p is equaled to k = 1. The suppressed power at lower k is because the in-homogenous 
of V p (r) at smaller k. The decay trend follows A; -5 / 3 and the fluctuations caused by the radial oscillatory function with random phases. 
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(a) 3D Perturbed Velocity 



(b) Kinetic Energy Power Spectrum of 3D perturbed field 



Figure 8. (a) 3D perturbed velocity field: The iso-surfaces show the magnitude of the perturbed velocity on physical scales of 10 12 
cm. (b) Corresponding energy power spectrum: The dotted red line shows the Kolmogorov spectrum, E(k) ~ k~ 5 / 3 . The peak of the 
Kolmogorov spectrum is adjusted to fit the data. The scale of H p is equaled to k = 1. Similar to 2D, The decay trend flows fc -5 / 3 and 
the fluctuations caused by the radial oscillatory function with random phases. 



speed, v is the fluid velocity and the timescale for burning is dtb, which is determined by both the energy 
generation rate and the rate of change of the abundances. 
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Figure 9. Total explosion energy as a function of resolution: the x-axis is the grid resolution and the y-axis is the total energy, defined 
to be the sum of the gravitational energy, the internal energy, and the kinetic energy. The total energy is converged when the resolved 
scale is close to 10 s cm. The right panel shows the zoom-in of the red box in left panel. 

5.1 Homographic Expansion 

As we have shown, grid resolutions of 10 8 cm are needed to fully resolve nuclear burning of our model. 
However, the star can have a radius of up to several 10 14 cm. This large dynamical range (10 6 ) makes it 
impractical to simulate the entire star at once while fully resolving all relevant physical processes. When 
the shock launched from the center of star, the shock traveling time scale is about few days which is much 
shorter than the Kelvin-Helmholtz time scale of the stars, about several million years. We can assume 
when the shock propagates inside the star, the stellar evolution of outer envelope is frozen. This allows us 
trace the shock propagation without considering the overall stellar evolution. Hence we instead begin our 
simulations with a coordinate mesh that encloses just the core of the star with zones that are fine enough 
to resolve explosive burning. We then halt the simulation as the SN shock approaches the grid boundaries, 
uniformly expand the simulation domain, and then restart the calculation. In each expansion we retain the 
same number of grids (see Figure [TO]) . Although the resolution decreases after each expansion, it does not 
affect the results at later times because burning is complete before the first expansion and emergent fluid 
instabilities are well resolved in later expansions. These uniform expansions are repeated until the fluid 
instabilities cease to evolve. There might be some possible sound waves generated from boundaries under 
such a setup. However, the normal SN shocks have a very higher mach number above 10 while traveling 
inside the star. The sound waves could not contaminate the burning/fluid instabilities domains before the 
shock reaches the boundary of simulation box. 

Most stellar explosion problems need to deal with a large dynamic scale such as the case discussed here. 
It is computationally inefficient to simulate the entire star with a sufficient resolution. Because the time 
scale of explosion is much shorter than the dynamic time of stars, we can follow the evolution of shock by 
starting from center of star and trace it until the shock breaks out the stellar surface. 



6 Conclusion 



Multidimensional stellar evolution and supernova simulations are numerically challenging because mul- 
tiple physical processes (hydrodynamics, gravity, burning) occur on many scales in space and time. For 
computational efficiency, ID stellar models are often used as initial conditions in 2D and 3D calculations. 
Mapping ID profiles onto multidimensional grids can introduce serious numerical artifacts, one of the most 
severe of which is the violation of conservation of physical quantities. We have developed a new mapping 
algorithm that guarantees that conserved quantities are preserved at any resolution and reproduces the 
most important features in the original profiles. Our method is practical for ID and 2D calculations, and we 
are now developing integral methods (an explicit integral approach instead of using volume subsampling) 
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Figure 10. Homographic expansion: In both panels, the yellow circle is the SN shock and the red region is the ejecta. The simulations 
begin with just the inner part of star and a higher resolution (left panel) for capturing fluid instabilities and burning. After the explosion 
occurs, we follow the shock until it reaches the boundary of the simulation box. We then expand the simulation domain, mapping the 
final state of the previous calculation onto the new mesh with a new ambient medium that is taken from the initial profile. 



that are numerically tractable in 3D. 

Multidimensional models give insight on fluid instabilities in supernova explosions that break the spher- 
ical symmetry of stars and mix their interiors. These instabilities originate from perturbations in the star 
prior to the explosion. Until now, these perturbations have been randomly seeded in 2D and 3D models 
with little or no physical basis. We present a new approach to seeding supernova models with physically 
realistic velocity perturbations like those found in the turbulent convective zones of massive stars. We 
find that the initial spectrum of the perturbations tends to be smeared out as they become nonlinear. Our 
approach can be applied to other multidimensional simulations of stellar explosions, especially those whose 
final outcomes are sensitive to the form of the initial perturbation; or the simulations of short duration, 
in which perturbations may not evolve fully nonlinear. 

Finally, we provide possible approaches to obtain the proper resolution for simulations that include 
both hydrodynamics and nuclear burning. Because the burning changes both the internal energy and 
composition of the fluid, we determine the physical scale for resolving burning with resolution tests and 
proper time steps by considering both hydro and burning. We apply a homographic expansion to bypass the 
numerical difficulties associated with the large range of dynamical scales in our problem. The algorithms we 
present can be applied to other multidimensional simulations besides stellar explosions in both astrophysics 
and cosmology. 
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